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ABSTRACT 

We study the evolution of a spherically symmetric density perturbation in the Modified 
Newtonian Dynamics (MOND) model applied to the net acceleration over Hubble flow. 
The background cosmological model is a A-dominated, low-fif, Friedmann model with 
no Cold Dark Matter. We include thermal processes and non-equilibrium chemical 
evolution of the collapsing gas. We find that under these assumptions the first low- 
mass objects (Af < 3 X IO^Mq) may collapse already for z ^ 30, which is in quite 
good agreement with the recent WMAP results. A lower value of gq would lead to 
much slower collapse of such objects. 
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1 INTRODUCTION 

Recent developments in cosmological observations have led 
to the so-called cosmological concordance model with Qi, 
around 0.03, Q,m (dark-|-baryonic) around 0.3 and JIa 
around 0.7. However, as the ACDM models are dominated 
by hypothetical vacuum energy and non-baryonic dark mat- 
ter contributions, some scientists are looking for different 
solutions. Perhaps the most intereresting alternative model 
is the M odified Newtonia n Dynamics model (MOND) pro- 
posed bv iMilgromI il983l) . It assumes that there is no non- 
baryonic dark matter (or it is negligible) and the apparent 
lack of matter is only due to the modification of dynamics or 
gravity for small accelerations (a ^ ao where ao is some con- 
stant). This model seems to work ver y well for spiral galax- 
ics and many other types of objects iSanders fc McGauebl 
L2OO2). However, it carries with it some unresolved problems 
(e.g. lack of covariance). 

The problem of the structure formation in MOND was 
also explored bv lSandcrs (200_1,) and Nusscr (2002), but they 
were interested in large cosmic structures and did not in- 
clude gas effects. It is a good approximation for large scales, 
but gas effects play a crucial role in the small-scale structure 
formation. 

IStachniewicz fc KutscheTal (I2OOIIJI studied possible im- 
plications of the MOND model for the formation of the very 
first objects in the Universe. However, there we had applied 
MOND to total gravitational acceleration, and this approach 
can lead to some paradoxes, e.g. time and the very occurence 
of collapse hardly depend on the initial overdensity. In this 
paper, we attempt a different approach; namely, we apply 
MOND to the net acceleration over the Hubble fiow only. 



2 MOND VS STANDARD THEORY OF 
LINEAR PERTURBATIONS 

If one wants to apply the MOND model to structure for- 
mation calculations, they encounter a number of difficulties. 
First of all, MOND is not a theory, but rather a phenomeno- 
logical model. In its present for m, MOND is incon sistent 
with General Relativity. Although iBekensteirJ ||2004) found 
a relativistic gravitation theory that gives MOND-like pre- 
dictions in the low-gravity limit, it needs to be verified. 

MOND is a model that modifies either dynamics or 
gravity (in this paper we assume the latter). It introduces 
a new fundamental scale, usually called ao. Gravitational 
fields much stronger than ao are identical to their Newto- 
nian limit qjv, and very weak fields are ^aoQN- According 



m an nmit qjy, ana very 
to lSanders fc VerheiierJ i 



19981) , the value of the fundamental 
acceleration scale is ao = 1.2 x 10~*cm/s^. More precisely, 
the strength of the gravitational field may be written as 

m( — )5'=ffiv> (1) 
Vao / 



where /i(a;) is some function that interpolates between these 
two extreme cases. This function is not specified in the 
model. We have decided to apply the function used by 
ISanders fc VerheiienI (Il998h : 



and, finally, 



9N 



(2) 



(3) 



where x — gisr/cLo- 
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The conseqences of MOND for c osmoloKy have not been 
studied in detail vet. lSandersI il998l) suggested that because 
in the early Universe the MOND radius is much lower than 
that of the horizon, the evolution of the scale factor is de- 
scribed by the s tandard Friedman n equations. Also in the- 
ory described bv lBekensteinI ||2004) cosmological models are 
similar to the Friedmann ones so we follow this assumption 
and study the formation of the first objects in the Universe 
with modified dynamics. 

We adopted this approach in lStachniewicz fc Kutscheral 
i2001lj) . but, as mentioned in the Introduction, this may 
lead to some paradoxes - in particular, even underdense re- 
gions of the Universe may collapse. In order to avoid that 
problem, we have tried another approach appliying MOND 
to the net acceleration/deceleration over the Hubble fiow. 
In this approach, only clouds with positive overdensity may 
stop their expansion and recoUapse, so the final formula for 
gravitational acceleration (or deceleration) is 



(4) 



where qh is the deceleration (or acceleration) due to the 
Hubble flow. 



2.1 Collapse of perfect gas 

If we take perfect gas instead of the pressureless fluid, the 
evolution will look different because the effects of pressure 
will moderate the recoUapse, especially for small systems. 
Since we assume spherical symmetry, Lagrangian coordi- 
nates are used. 

In the Newtonian case, dynamics is governed by the 
following equations: 



dM 




(5) 


dr 


— 4nr^g, 


dr 




(6) 


lit 


= V, 


dv 


_ 4^,2 dp GM{r) 


(7) 


It 


dM r2 ' 


du 


p dg ^ A 


(8) 


IE 


g^ dt g' 



where r is the radius of a sphere of mass M , it is the internal 
energy per unit mass, p is pressure and g is mass density. 
Here, Eq.l|KJ is the continuity equation, @ and lO give ac- 
celeration and accounts for energy conservation. The last 
term in Eq.ljHl describes gas cooling/heating, with A being 
energy a bsorption (emission) rate per unit volume, given in 
detail in lStachniewicz fc Kutscheral ||2Q^3)- 
We use the equation of state of perfect gas 



(7 - l)gu, 



(9) 



where 7 = 5/3, as the primordial baryonic matter af- 
ter recombination is assumed to be composed mainly of 
monoatomic hydrogen and helium, with the fraction of 
molecular hydrogen H2 always less than 10~^. 

In the case of modified gravity, equation Q will look 
somewhat different: 



dv 
IE 



dM 
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aor^ 



(10) 



where f(x) is inverse to function ^i{x) mentioned before, 
asymptotically equal to x for a: 2> fflo and to y/oox for x <C 
ao, while qh may be expressed as 



(11) 



Here, Ho is the current value of the Hubble parameter, fit,, 
fir and f^A are the current fractions of baryons, radiation 
and dark energy in terms of the critical density of the Uni- 
verse, and z is the redshift. 

2.2 Chemical reactions and thermal effects 

In our calculations, we include all of the relevant thermal 
and chemical processes in the primordial gas. In this pa- 
per, we have taken into account nine species: H, H~, H^, 
He, He^, He^""", H2, Hj and e~. The abundance of various 
species changes with time due to chemical reactions, ioniza- 
tion and dissociation photoprocesses. The chemical reactions 
include such processes as ionization of hydrogen and helium 
by electrons, recombination of ions with electrons, formation 
of negative hydrogen ions, formation of H2 molecules, etc. A 
full list of the relevant chemical reactions and appriopriate 
formulae is given in ,, Stachniewic z fc Kuts chera (2001a). 

The time evolution of the number density of component 
Ui is described by the kinetic equation: 



dui _ \ - \ - 
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(12) 



The first component on the right-hand side of this equation 
describes the chemical reactions, and the other one accounts 
for photoionization and photodissociation processes. Coeffi- 
cients kim are reaction rates, quantities k„ are photoioniza- 
tion or photodissociation rates, and aimi and bji are numbers 
equal to 0, ±1 or ±2 depending on the reaction. All reaction 
rates, as wel l as photoionization and photodisso ciation rates 
are given in IStachniewicz fc Kutschem (12001311 . 

The cooling (heating) function Acooi includes effects 
of collisional ionization of H, He and He^, recombina- 
tion to H, He and He"*", collisional excitation of H and 
He+, Bremsstrahlung, Compton cooling and cooling by de- 
excitation of H2 molecules. The formulae for the heat- 
ing/cooling contributions o f vario us processes are given in 
IStachniewicz fc KutscheTal (l2001afl . 



3 CODE USED IN THE SIMULATIONS AND 
INITIAL CONDITIONS 

In the simulations we used t he code described in 
IStachnie wiczj£jCutschgra|_j|20 l a|) . b ased on those pre- 
sented bv~ lThouO^Weinbe^ ~ ] 1995^ and iHaiman et alJ 

(1996). This is a standard, one-dimensional, second-order 
accurate, Lagrangian finite-difference scheme. The only 
changes consist in the modification of gravity, with the dark 
matter fraction Qdm set to zero. However, it was necessary 
to make significant changes in the initial conditions. 

First of all, we started our calculations at the end of the 
radiation-dominated era. For fit = flm = 0.02/ h^, Zea = 
485 as given by the formula provided by iHu fc EisensteinI 
(1998), Zeq = 2.50 X lO^f^o/i^e^T, where 82.7 = T^/2.7K, 
assuming h = 0.72 and r^=2.7277 K. We assumed that, like 
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Figure 1. Matter power spectrum. 
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Figure 2. Turnaround times for shells of various masses and 
initial overdensities 



in the standard cosmology, initial overdensities may grow 
only in the matter-dominated era. We used our own code 
to calculate the initial chemical composition and initial gas 
temperature. When compared with the results obtained by 
iGalU fc Pall3 il998 l). ours are very similar - the agreement 
is at the level of 10-20%. The difference is probably due to 
the fact that Galli and Palla had included more species (e.g. 
deuterium and lithium), and some reaction rates they used 
had been somewhat different. 

It is necessary to specify the initial baryonic matter 
perturbations, 5. We have decided to take values of the or- 
der of the contribution to the baryonic matter power spec- 
trum at a given scale and redshift. The problem is that 
at these scales the power spectrum is known very roughly, 
as there are no observational data and the relevant phys- 
ical processes are not known very well. The power spec- 
trum may be calcula t ed usi ng the CM BFAST program by 
ISeliak fc Zaldarriaeal Jl996D A{k) = y/'lndlorrnkHf{k,z)^, 
where tf is the transfer function, k is the wavenumber 
and the value of normalization constant d^^^^ is obtained 
from CMBFAST. We also included the output of the CMB 
anisotropy program written by N. Sugiyama, received from 
the author. The results are shown in Fig. 

The most interesting for us are the comoving scales 
around 10 kpc. For these scales, both outputs differ signif- 
icantly by a factor of 100. We decided to take our initial 
overdensity from in between, i.e. equal to 10~^. As we will 
show, these differences do not affect our results significantly. 

We apply the initial density profiles in the form of a 
single spherical Fourier mode, also used by iHaiman et alJ 



Qbir) = ilbQcil + S 



sin kr ^ 
kr 



(13) 



where Qc is the critical density of the Universe, Qc = 
SHVSttG, with H being the actual value of the Hubble pa- 
rameter. 

For this profile, there exist two distinguished radius val- 
ues, namely Rq and Rz, which correspond to the first zero 
and the first minimum of the function sin(fcr)/fer, respec- 
tively. Inside the sphere of radius Ro = Tv/k which contains 
mass Mo, the local density contrast is positive. The mass Mo 
and the radius _Ro will be referred to as the cloud mass and 



the cloud radius, respectively. The local density contrast is 
negative for R^ > r > Rq, with average density contrast 
vanishing for the sphere of radius R^ = 4.49341/fc and mass 
Mz . According to the gravitational instability theory in the 
expanding Universe, the shell of radius R^ will expand to- 
gether with the Hubble fiow, not undergoing any additional 
deceleration. This is why we regard this profile as very con- 
venient in numerical simulations. Namely, it eliminates the 
numerical edge effects and the mentioned shell simply fol- 
lows the Hubble expansion of the Universe. Thus, it can be 
regarded as the perturbation boundary, with mass A'lz re- 
ferred to as the bound mass. 

It is worth noting that for radii not greater than 3/4_Ro 
this profile is very similar to the Gaussian profile 



gi{r) = QiQc 



1 + 5i exp 



2Rj 



(14) 



with Rf = l/2Ro. 

As the initial velocity, we use the Hubble velocity: 

v{r) = Hr. (15) 



4 PRESSURELESS COLLAPSE IN MOND 

We decided to describe the case of a pressureless collapse 
in MOND as this makes it easier to understand why the 
results shown in the next section hardly depend on the initial 
overdensities. It is much easier to trace shells since each 
shell may be treated separately, so it is enough to trace the 
evolution of a single shell with some initial overdensity and 
constant mass inside. Equation 1101 is reduced to 



dv ^ ( GM(t 

— = -9h - aof {- 



dt 



aor^ 



iE. 

ao 



(16) 



where the expression for gn was shown in Eq. 1111 We wrote 
a small program calculating the turnaround times for shells 
of various masses and initial overdensities. It starts at the 
beginning of the matter-dominated era with the initial ve- 
locity being the Hubble velocity (v = Hr). The results are 
shown in Fig. |5| 
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As we can see, for initial overdensities of 10~* and lower, 
the results are almost indistinguishable and not very differ- 
ent from those for 10~*. This looks very odd, but it reflects 
the fact that a smaller overdensity means that the shell is 
very deep in the MOND regime and its behaviour is strongly 
nonlinear. Now let us assume that a shell is always in the 
MOND regime. For pure MOND Eq.[Tnitakes the form 



dv 
IE 



-9H 



GM{r) 



9H 



(17) 



It is a nonlinear equation, especially in the low over- 
density limit. This becomes obvious if we change variables. 
Instead of r, let us put {R — r), where R is the radius corre- 
sponding to the mass M in a homogenous Universe. In other 
words, r is related to the actual overdensity 5 and if it is 
much lower than R, it may be expressed as 5 = 3r/_R. Given 
that, equation 1171 mav be written as 



/ M r 



(18) 



Now we may express mass in terms of the cosmological pa- 
rameters, redshift and R, and we get 



Redshift z is a function of time, and if we 
dr = [|r2f,aoHo (z -|- lf]^''^dt, we finally obtain 



dfr_ 
dr2 



(19) 
put 

(20) 



5 RESULTS 

We performed seventeen runs. The flrst eight were for the 
'standard' value of ao (1.2 x 10~*cm/s^), various masses of 
the cloud (10^ Mq, 3 x 10^ Mq, IO^Mq and 3 x lO^M©) and 
initial overdensities (10~^ and 10~*, the next three were for 
some lower value of ao (1.2 x 10~^cm/s^), initial overdensity 
10"^ and masses of the cloud equal to 3 x 10^ Mq, 10'' Mq 
and 3 x 10* Af©, while the last one was for the 'standard' 
value of ao but with no H2 cooling. The results are shown 
in Figs. OlTTl 

Most of these flgures show trajectories of shells enclos- 
ing 7%, 17%, 27% ... 97% of the total mass; Figs.EHfor the 
'standard' ao and 10~^ overdensity and Figs. I9I11I for the 
'low' ao, 10~^ overdensity, respectively. In addition, Fig. Q 
shows temperatures of the same shells as in Fig. I^and Fig. 
|S| shows the chemical evolution for shells enclosing 0.12Mz, 
both figures showing runs for the 'standard' ao, 10~® over- 
density and M = 10'' M© (the results for the other masses 
are very similar). The results for 10~* are almost idistin- 
guishable and the run with greater mass {M = IQ^Mq) but 
no H2 cooling showed no collapse. 

The results show that: 

• as expected, the difference in behaviour between clouds 
with 10~^ and 10~* overdensities is very tiny, thus the re- 
sults are much less sensitive to the initial density contrast 
than in CDM 

• more massive clouds collapse faster 

• speed of collapse depends very strongly on ao 
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Figure 3. Shell trajectories for the 'standard' ao, 10 ^ overden- 
sity, M = lO^A/©. 



0.1 



J; 0.01 

DC 



0.001 




100 



10 



z+1 



Figure 4. Shell trajectories for the 'standard' ao, 10 ^ overden- 
sity, M = 3 X 10^ M0. 
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Figure 5. Shell trajectories for the 'standard' ao, 10 ^ overden- 
sity, M = 1O*M0. 
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Figure 6. Shell trajectories for the 'standard' ao, 10 ^ overden- 
sity, M = 3 X IO^Mq. 
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Figure 9. Shell trajectories for the 'low' ao, 10 ^ overdensity, 
M = 3 X IO^Mq. 
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Figure 7. Shell temperaturs for the 'standard' ao, 10 ® overden- 
sity, M = IQ^Mq. 
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Figure 10. Shell trajectories for the 'low' ao, 10 ® overdensity, 
M = 10* Mq. 
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Figure 8. Chemical evolution for the 'standard' ao, 10 ^ over- 
density, M = IO^Mq. 
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Figure 11. Shell trajectories for the 'low' ao, 10 ® overdensity, 
M = 3 X IO^'Mq. 
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• like in the CDM models, H2 cooling is necessary to col- 
lapse 

• like in the CDM models, due to the adiabatic cool- 
ing/heating, the shell temperature behaviour is opposite to 
that of the shell radii; then H2 cooling becomes important 
and the shells collapse; a higher mass of the cloud means 
higher virial temperature and faster collapse 

• chemical evolution is quite typical but because the col- 
lapse is very violent, it speeds up chemical reactions; how- 
ever, the final abundances of various species are not very 
different from the predictions of the other models, e.g. the 
final abundance of H2 is of the order of 10~^. 

It may appear that during the collapse, shell tempera- 
tures do not fall as deep as in CDM mod els, but this is a nu- 
merica l artefact only. The reason is, as in lThoul fc Weinberg! 
il995l) . that in our code - if a shell falls below some 'small' 
radius Vc - we treat it as 'collapsed': we artificially stop 
its dynamical and chemical evolution, assigning some small 
radius, freezing the temperature and including its gravita- 
tional field only. A collapse in MOND is very violent, so gas 
shells have not enough time to cool down until their tem- 
perature is frozen by the code. 

We performed many other runs as well, e.g. with initial 
overdensities greater than 10"*, but the results were very 
similar. The reason for this is that, as shown in the previous 
section, in the absence of gas pressure, any positive initial 
overdensity lower than some reasonable value leads to a very 
similar turnaround time. Because at this stage of evolution 
the gas effects are not yet very important, this also applies 
to the evolution with all physical processes included. This 
means that in the approximation the turnaround time de- 
pends on mass. Then the cloud starts its initial collapse, but 
it gets virialised and the time between the virialisation and 
the final collapse depends on mass again, i.e. the bigger the 
mass, the faster the collapse. 

The Wilkinson MAP results suggest that reioniza- 
tion occured at around z ^ 20. This means that the 
first bound objects must have been formed even ear- 
lier, perhaps at around z ^ 30 . Our recent simulations 
iStachniewicz fc Kutscheral I2OO3I) show that if we assume 
the ACDM models and take the recent estimates of Qm 
and JIa, a direct formation of low-mass objects that could 
possibly reionize the Universe before 2 ~ 10 is very un- 
likely. Moreover, we doubt whether the inclusion of possible 
fragmentation of greater clouds could speed up the collapse 
enough - even if some low-mass cloud has greater overden- 
sity than the directly forming ones, it still needs some time 
to cool down. 

In contrast, MOND seems to provide a good way to 
solve that problem. For the 'standard' value of ao clouds of 
mass 3 x 10^ Mq or heavier may collapse at around z ~ 30, 
so that they or their cores may have formed the first stars 
and quasars. For lower ao, only objects of mass 3 x 10'' M0 or 
greater may be formed directly before z ~ 30. We think this 
favours the 'standard' value but one would need to perform 
full 3-D simulations to give a more definite answer. 

It may be surprising that in our simulations more mas- 
sive objects form earlier, but it is due to the gas pressure, 
which is more important for less massive objects. If we skip 
that term in Eq.|7| less massive objects stop their expansion 
and recoUapse faster than the more massive ones. However, 



this is not the case for the Large Scale Structure: on these 
scales, gas effects may be neglected, and although in general 
the MOND effects are smaller, the rms fluctuations decrease 
so that such objects stop their expansion more slowly. 



6 CONCLUSIONS 

If our assumptions about MOND are correct, its predic- 
tions seem to be more consistent with t he the early reioniza- 
tion suggested by the WMAP results jBennett et alJl2003l: 
ISpergel et al.ll20o3) than the ones of the most recent ACDM 
models. This does not prove that MOND is correct and 
ACDM is not. However, this suggests that cosmologists 
should perhaps pay more attention to MOND, because it 
seems to be an interesting alternative to models with non- 
baryonic dark matter. 
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